#include <vector>
#include <random>
#include <algorithm>


//------------------------------------------------------------------------------
// Given a table of numeric categories (e.g., ages) and a table of frequencies
//   draws a random category (age).
//------------------------------------------------------------------------------
class CInvDistribution {

 private:

   std::vector<double> V;
   std::vector<double> Categories;

   std::uniform_real_distribution<double> dist;

   std::mt19937 *generator;

   double Q;

   //---------------------------------------------------------
    int find_idx (double freq) {

        if (V.size() < 1) return -1;

        int k = 0;
        for (k = 0; k < V.size(); ++k) {
           if (V[k] > freq) break;
        }

        return k;
    }

 public:
 
   CInvDistribution (std::mt19937 *p) : dist(0.0,1.0) {
      generator = p;
   }

   CInvDistribution (std::vector<double> p) : dist(0.0,1.0) {
      prepare(p);
   }

   CInvDistribution (std::vector<double> p,
                     std::vector<double> c) : dist(0.0,1.0) {

      if (p.size() != c.size())
         return;

      prepare (p);

      Categories = c;
    }


   void prepare (std::vector<double> v) {

      Q = 0.0;

      V.clear ();

      for (auto a:v) {

         Q = Q + a;

         V.push_back (Q);

      }
   }

   //---------------------------------------------------------
   void prepare (std::vector<double> v, std::vector<double> c) {

      if (v.size () != c.size())
         return;

      prepare (v);
      Categories = c;
   }

   //---------------------------------------------------------
   int idraw (void) {

      return find_idx ( Q * dist(*generator) );
   }

   //---------------------------------------------------------
   double ddraw (void) {

      double a;
      double r = Q * dist(*generator);

      int k = find_idx ( r );

      if (k == 0) {
         
         a = 1.0 -  ((V[0] - r) / V[0]); // assumindo idade mínima é zero (não negativa)

         return Categories[1] * a;
      }

      //if (k >= V.size() - 1) return Categories[Categories.size()-1]; // saturação?

      a = (r - V[k-1]) / (V[k] - V[k-1]);

      return Categories[k-1] + a * (Categories[k] - Categories[k-1]);

    }
};


//------------------------------------------------------------------------------
// Given a table of numeric categories (e.g., ages) and a table of frequencies
//   returns a probability given a category.
//------------------------------------------------------------------------------

class CProbTable {

  private:

   std::vector<double> V;
   std::vector<double> Categories;

  public:

   CProbTable () {;};

   //---------------------------------------------------------
   void prepare (std::vector<double> v, std::vector<double> c) {

      if (v.size () != c.size())
         return;

      V          = v;
      Categories = c;
   }

   //---------------------------------------------------------
   double Probability_Of (double cat) {

      if (Categories.size() < 1) return -1;

      int k = 0;
      for (k = 0; k < Categories.size(); ++k) {
         if (cat < Categories[k]) break;
      }

      if (k == 0) return V[0];

      if (cat >= Categories[Categories.size()-1]) return V[V.size()-1];

      double dx = Categories[k] - Categories[k-1];
      double dy = V[k]          - V[k-1];

      double alpha = dy / dx;

      return V[k-1] + alpha * (cat - Categories[k-1]);

   }
};
